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Abstract. In this article algorithmic methods are presented that have essen- 
tially been introduced into computer algebra systems like Maple or Mathematica 
within the last decade. The main ideas are due to Stanley and Zeilberger. Some 
of them had already been discovered in the last century by Beke, but because 
of their complexity the underlying algorithms have fallen into oblivion. We give 
a survey of these techniques, show how they can be used to identify transcen- 
dental functions, and present implementations of these algorithms in computer 
algebra systems. 

1 Algebraic Representation of Transcendental 
Functions 

How can transcendental function be represented by algebraic means? To give 
this question another flavor: What is the main difference between the exponen- 
tial function f(x) = e x and the function g{x) = e x + \x\/10 1000 , that makes / 
an elementary function, but not <?, although / and g are numerically quite close 
on a part of the real axis? 

Or let's consider an example of discrete mathematics: Why is the factorial 
function a n — n\ considered to be the most important discrete function, and 
not b n = n! + n/10 1000 or any other discrete function? 

Although these examples refer to the most important continuous and discrete 
transcendental functions, oddly enough the answers to the above questions are 
purely algebraic: The exponential function / is characterized by any of the 
following algebraic properties: 

1. / is continuous, /(l) = e, and for all x,y we have f(x + y) = f{x) ■ f(y); 

2. / is differentiable, f'(x) = f(x) and /(0) = 1; 

oo 

3. / 6 C°°, f(x) — Yl a n x n with ao = 1, and for all n > we have 

71=0 

(n + 1) On+i = a„; 

and the factorial function a n is represented by any of the following algebraic 
properties: 

4. ao = 1, and for all n > we have a„_|_i = (n + 1) a n ; 

oo 

5. the generating function f(x)= ^2 a n z n satisfies the differential equation 
x 2 f'(x) + (x- l)f(x) + 1 = with the initial condition /(0) = 1. 
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(Note here that one could argue that property (1.) is not algebraic since the 
symbol e is needed in the representation.) I do not know any method to represent 
transcendental functions using functional equations, such as property (1.), but 
I will show, why and how the other properties can be suitable for this purpose, 
being mainly concerned with properties (2.) and (4.). In § [| we consider, how 
these representations can be viewed as purely polynomial cases. 

Observe that the "generating function" of the factorial function is conver- 
gent only at the origin, and therefore must be considered as a formal series. In 
particular, a "closed representation" (whatever that should mean) of the gen- 
erating function cannot be given. But this is not the main issue here. Rather 
than working with the generating function itself, it is much better to work with 
its differential equation which is purely algebraic (in fact, it is purely polyno- 
mial). The same argument applies to the exponential and factorial functions 
themselves. Rather than working with these transcendental objects, one should 
represent them by their corresponding differential and recurrence equations. 

The given properties are structural statements about the corresponding func- 
tions. Any small modification (even changing the value at a single point) de- 
stroys this structure. For example, the function g(x) = e x + |a;|/10 1000 cannot 
be characterized by a rule analogous to one of the properties (2.)-(3.). On the 
other hand, the function h(x) = e x +a;/10 1000 can be represented by the differ- 
ential equation (x— 1) h"(x) — x h'(x) +h(x) = with the initial values h(0) = 1 
and h'(0) = 1 + 10- 1000 . 

Therefore, the special (and common) fact about the exponential and fac- 
torial functions is that they both satisfy a differential or recurrence equation, 
respectively, that is homogeneous, linear, of order one, and has polynomial co- 
efficients. 

We can generalize this observation |42| : A continuous function of one variable 
f{x) is holonomic, if it satisfies a homogeneous linear differential equation with 
polynomial coefficients; we call such a differential equation also holonomic. 

By linear algebra arguments, Stanley |3r| showed that sums and prod- 
ucts of holonomic functions and the composition with algebraic functions also 
form holonomic functions. This can be seen as follows: Assume / and g sat- 
isfy holonomic differential equations of order n and m, respectively. We con- 
sider the linear space Lf of functions with rational coefficients generated by 
/, /', /", . . . , f( k ', . . .. Since /,/',.-■, f^ n ' are linearly dependent by the given 
holonomic differential equation and since by differentiation the same conclusion 
follows for /', /", . . . , and so on inductively, the dimension of Lf is < n. 

Similarly L g has dimension < m. We now build the sum Lf + L g which is of 
dimension < n + m. As f+g, (/+ g)' , ■ ■ ■ , (/+<?) ^ , • • • are elements of Lf + L g , 
arbitrary n + m + 1 many of them are linearly dependent. In particular, f + g 
satisfies a holonomic differential equation of order <n + m. 

Similarly the product and composition cases can be handled. Note that 
the above proof provides a construction of the resulting holonomic equation 
by linear algebra techniques. It is remarkable that 100 years ago, Beke 
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already described these algorithms to generate holonomic differential equations 
for the sum and product of / and g from the holonomic differential equations 
of / and g. Hence, he had discovered algorithmic versions of Stanley's results! 

Analogously, a discrete function (sequence) of one variable is called holo- 
nomic, if it satisfies a homogeneous linear recurrence equation with polynomial 
coefficients. Such a recurrence equation is also called holonomic. Sums and 
products of discrete holonomic functions are again holonomic, and there are 
similar algorithms to calculate representing holonomic recurrence equations (s. 

©, ©)■ . 

A function 

oo 

71=0 

represented by a power series is holonomic if and only if the corresponding power 
series coefficient a n is a holonomic sequence. The holonomic equations for f{x) 
and a n can be converted equating coefficients. 

Note that these algorithms were implemented by Salvy and Zimmcrmann in 
the gfun package of Maple's share library j34|. 1 wrote a Mathematica imple- 
mentation, SpecialFunctions, to be obtained by World Wide Web from the ad- 



dress f tp : //f tp . zib-berlin. de/pub/UserHome/Koepf /SpecialFunctions. 



Examples of this implementation will be given later. 



2 Identification of Transcendental Functions 

Note that the notion of holonomy provides a normal form for a suitably large 
number of transcendental functions, which can then be utilized for identification 
purposes. The holonomic equation of lowest order corresponding to a holonomic 
function constitutes such a normal form. Once we have calculated the normal 
form of a holonomic function, the latter is identified: Two holonomic functions 
are identical if and only if they have the same normal form, and satisfy the same 
initial conditions. 

But also without having access to the lowest order holonomic equations, 
one can check whether two holonomic functions agree, since (by linear algebra, 
e.g.,) it is easy to see whether two holonomic equations are compatible with 
each other. 

Therefore, we may ignore that e x ,sin:z:,cosx,arctan2;,arcsin:r and others 
form transcendental functions, and take only their holonomic differential equa- 
tions /' = /, /" = -/, /" = -/, (1 + x 2 )f" + 2xf = 0, (x 2 - l)f" + xf = 
etc. into account. From these differential equations, corresponding differen- 
tial equations for sums and products can be generated by the above mentioned 
technique, using only polynomial arithmetic and linear algebra. For example, 
the function f(x) = arcsin 2 x yields (x 2 - 1)/"' + 3xf" + f = 0. Note, how- 
ever, that in the given case one can get even more: The resulting holonomic 
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differential equation is directly equivalent to the holonomic recurrence equa- 
tion n(l + n)(2 + n)a n +2 = n 3 a n for the coefficients a n of the Taylor series of 

oc 

arcsin 2 x = ^] a n x n , and since this holonomic recurrence equation fortunately 

contains only the two terms a n +2 and a n , it can be solved explicitly, and leads 
to the representation 

4™n' 2 

arcsin 2 x — > -. -— -x 2n+2 

^ (l + n)(l + 2n)! 

(compare @, @, 

Note that not only a function like the Airy function Ai (x) (s. e. (10.4)) 
falls under the category of holonomic functions, since it satisfies the simple 
holonomic differential equation /" — xf = 0, moreover the classical families of 
orthogonal polynomials^ and many other special functions form holonomic func- 
tions Q. These depend on several variables, and we will discuss this situation 
in § | 

On the other hand, there are functions that are not holonomic, like the 
tangent function tana: (s. |36j, |25[]). The identification problem for expressions 
involving nonholonomic functions can only be treated after preprocessing the 
input. If, for example, we want to verify the addition formula for the tangent 
function 

tan x + tan y 

tan (x + y) = 

1 — tan x tan y 

by the given method, then we have to replace all occurrences of the tangent 
function by sines and cosines (which are holonomic) using the rewrite rule 
tan a; = sin xj cos x . We can then generate a polynomial equation by multi- 
plying both sides by the common denominator. This procedure results in the 
equivalent representation 

(cos a; cosy — sin x sin y) sin (x + y) = (cosy sin x + cos x sin y) cos (x + y) (1) 

which is easily proved since the algorithms generate the common holonomic 
differential equation f"(x) + 4/'(x) = with respect to x (or the common 
holonomic differential equation f"(y) + 4/'(y) = with respect to y) for both 
sides of (0) where the common initial values are /(0) = cosy siny, and /'(0) = 
cosy 2 — siny 2 . Assume that for the initial value functions we had obtained 
different representations (e.g. cosy siny and sin(2y)/2). These could be verified 
by the same technique. 

In the Mathematica package SpecialFunctions (s. also |^2|), the procedure 
HolonomicDE [f ,x] calculates the holonomic differential equation of / with re- 
spect to the variable x using the known holonomic differential equations of the 
primitive functions, and the sum and product algorithms by recursive decent 



x As families of orthogonal polynomials they are not polynomials! 
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through the expression tree. Here we call a function primitive if it is rational, 
or whenever we use a separate symbol for it and a holonomic differential equa- 
tion is known. Therefore the above mentioned functions (besides the tangent 
function) are primitive. 

The examples given are governed by the following Mathematica session: 

In[l]:= «SpecialFunctions' 

In [2] := HolonomicDE[ArcSin[x] ~2,x] 

(3) 

Out [2]= F'[x] + 3 x F"[x] + (-1 + x) (1 + x) F [x] == 
In[3]:= DEtoRE[ , /.,F,x,a,n] 
3 

Out [3]= n a[n] - n (1 + n) (2 + n) a[2 + n] ==0 

In [4] := Series [ArcSin[x] ~2,{x,0}] 

k 2 + 2 k 2 
4 x k! 

Out [4]= Sum[ , {k, 0, Infinity}] 

(1 + k) (1 + 2 k) ! 

In [5] := HolonomicDE [AiryAi [x] ,x] 

0ut[5]= -(x F[x]) + F"[x] == 

In [6] := HolonomicDE [AiryAi [x] ~2, x] 

(3) 

0ut[6]= 2 F[x] + 4 x F'[x] - F [x] == 

In [7] := HolonomicDE [Sin [x+y] * (Sin [x] Sin [y] -Cos [x] Cos [y] ) ,x] 
(3) 

Out [7]= 4 F> [x] + F [x] == 

In [8] := HolonomicDE [Cos [x+y]* (Sin [x] Cos [y] +Cos [x] Sin [y] ) ,x] 
(3) 

Out [8]= 4 F' [x] + F [x] == 
In [9] := HolonomicDE [Cos [y] *Sin[y] ,y] 
(3) 

Out [9]= 4 F> [y] + F [y] == 



In [10] := HolonomicDE[Sin[2y]/2,y] 
Out [10]= 4 F[y] + F"[y] == 

One difficulty that may arise with the method described is that in some in- 
stances the sum and product algorithms will not generate the holonomic differ- 
ential equation of lowest order, as in the above example for cosy sin y. In this 
case, the normal form property is lost. In fact, the sum algorithm calculates 
a holonomic equation that is valid for any linear combination af + bg rather 
than the particular given sum / + g. As a simple example, we consider the sum 
y/1 + x + ^l +x satisfying the first order differential equation 

2 (2 + x) (1 + x) F'(x) - xF{x) = . 

This differential equation can be found using a method given in [^0| |2l| , whereas 
the sum algorithm generates the second order differential equation 

4 (1 + xf F"(x) + 4 (1 + x) F'{x) - F(x) = . 

The reason for the existence of a differential equation of lower order is due to 
the fact that the ratio of the two summands \/l + x and ,} forms a rational 
function. 

Similarly, the sum of two consecutive Legendre polynomials P n {x) + P n +\(x) 
satisfies the second order differential equation 

(x -l)(x+ 1) F"{x) + {x + 1) F'(x) - (n + l) 2 F(x) = , 

whereas the sum algorithm generates the differential equation 

= {x-lf(l + x) 2 F""(x) + 8(x-l)x(l + x) F"'(x) 
+2 (-2 + 2 n + n 2 + 6 x 2 - 2 n x 2 - n 2 x 2 ) F"(x) 
-4 n (2 + n) x F'(x) + n (1 + n) 2 (2 + n) F(x) 

of fourth order, which is also valid for the difference P n {x) — P n+ i(x) and for 
any other linear combination. 

For the verification of identities, this is not an important issue, since the 
compatibility of two holonomic equations can be easily checked. This situation 
is similar to proving a rational identity by pure polynomial arithmetic without 
gcd computations (after having multiplied through by all denominators), and is 
actually equivalent to a noncommutative polynomial division, see § |]. 

In the case that the normal form is needed for a particular problem, a fac- 
torization algorithm can be used, s. § [| 

For the discrete functions, the situation is quite similar. We call a function 
primitive whenever we use a separate symbol for it and a holonomic recurrence 



G 



equation is known. To these primitive functions, we add the rational functions 
and the functions 



(mn + b)\, — T (meQ), and a n (2) 

[mn + by. 

whose holonomic recurrence equations are known, as primitive functions with 
respect to the variable n. We consider the factorial function to be equivalent 
to the r function T (a + 1) — al, and declare binomial coefficients etc. also via 
factorials. From the holonomic representations of the primitive functions the 
holonomic equations for all sums and products can be established. E. g. the two 
equations 

(n-k + l) 2 F(n + 1, k) - (1 + n) 2 F(n, k) = (3) 

and 

(k + l) 2 F(n, fc + 1) - (n - k) 2 F(n, k) = (4) 

for F(n, k) — (^) • Whereas these are simple consequences of the representation 
of F (n, fc) by factorials, the given procedure can be applied, for example, to the 
more complicated function F(n, fc) = " !H fc fc! to generate the two holonomic 
equations 

nF(n + 2, fc) - (1 + 3n + n 2 )F{n + 1, fc) + (1 + n) 2 F(n, fc) = 

and 

k(2+k) 2 F(n, fc+2)-(l+fc)(l+3fc+fc 2 )(3+3fc+fc 2 )F(n, fc+l)+fc(l+fc) 3 F(n, fc) = . 

Note that the given approach also covers all kinds of orthogonal polynomials 
and special functions with respect to their discrete variables, see § ^. 

In our Mathematica implementation SpecialFunctions, the procedure 
HolonomicRE [a,n] calculates the holonomic recurrence equation of a n with re- 
spect to the variable n taking the known holonomic recurrence equations of the 
primitive functions into account, and using the sum and product algorithms by 
recursive decent through the expression tree. The above examples are governed 
by the following Mathematica session: 

In[ll] := HolonomicRE [Binomial [n,k] "2, n] 

2 2 
Out [11]= (1 + n) a[n] - (1 - k + n) a[l + n] == 

In[12] := HolonomicRE [Binomial [n,k] ~2,k] 

2 2 
Out [12]= (-k + n) a[k] - (1 + k) a[l + k] ==0 
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In[13] := HolonomicRE[(n!+k! -2)/k,n] 

2 2 
Out [13]= (1 + n) a[n] + (-1 - 3 n - n ) a[l + n] + n a[2 + n] == 

In [14] := HolonomicRE[(n!+k! ~2)/k,k] 

3 

Out [14]= k (1 + k) (3 + k) a[k] - 

2 2 

> (1 + k) (1 + 3 k + k ) (3 + 3 k + k ) a[l + k] + 

2 

> k (2 + k) a [2 + k] == 



3 Hypergeometric Sums 

Rather than having functions given as finite sums and products of primitive 
expressions, an important class of functions, particularly in analysis and com- 
binatorics, is given by infinite sums of products of terms of the form (|2|) 

«(n) = X>(n,fc). (5) 

Then F[n, k) is an (m, i)-fold hypergeometric term. That is, both F(n + 
m, fc) / F(n,k) and F(n,k + I) / F(n,k) are rational functions with respect to 
n and k for a certain pair (m, I) € N 2 . For example, by this is valid for 

F(n, k) = (^) with m = I = 1. We assume moreover that the sums (|^) have 
finite support, i.e., they are finite sums for each particular n G N. 

A modification |23| of the (fast) Zeilberger algorithm see also f27f , 

and returns a holonomic recurrence equation valid for s(n). Zeilberger's 
algorithm is based on a decision procedure for indefinite summation due to 
Gosper jl7j. In our example case, Zeilberger's algorithm finds the holonomic 
recurrence equation (1 + n) s(n + 1) = 2(1 + 2n) s(n) for s(n) = (fe) 2 = 

fcGZ 

^2 (fc) which fortunately has only two terms. Therefore, we are led to the 

fe=0 

representation 




Even though, in general, the resulting recurrence equation has more than two 
terms, this holonomic equation contains very important structural information 
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about s(n). This may be used to show that a certain family of polynomials is 
orthogonal or not Q , and can be an interesting property for numerical purposes 
(compare Q). 

In particular, as described in the last section, the generated structural in- 
formation can be used for the identification of a transcendental function that 
is given as sum (@) . Note that sums of type (0) in general form transcendental 
functions with respect to the discrete variable n. 

For example, to check the identity (compare |f37f ) 

fc=0 v 7 fc=o v 7 v 7 

which is nontrivial since for n = 1 it reads 1 + 1 = + 2, we need only to show 
that both sums s(n) satisfy the common recurrence equation 

(n + 2) 2 s(n + 2) - (16 + 21n + 7n 2 )s{n + 1) - (n + l) 2 s(n) = (7) 

which is the result given by Zeilberger's algorithm. We also have the same initial 
values s(0) = 1 and s(l) = 2, so we are done. 

In Mathematica these computations are done by 

In [15] := HolonomicRE [Sum [Binomial [n,k] ~2, {k,0,n}] ,n] 

0ut[15]= -2 (1 + 2 n) a[n] + (1 + n) a[l + n] ==0 

In [16] := HolonomicRE [Sum [Binomial [n,k] "3, {k,0,n}] ,n] 

2 2 
Out [16]= -8 (1 + n) a[n] + (-16 - 21 n - 7 n ) a [1 + n] + 

2 

> (2 + n) a[2 + n] == 

In[17] := HolonomicRE [Sum [Binomial [n,k] ~2*Binomial [2k, n] ,{k,0,n>] ,n] 

2 2 
Out [17]= -8 (1 + n) a[n] + (-16 - 21 n - 7 n ) a [1 + n] + 

2 

> (2 + n) a[2 + n] == 

Note that the example shows that transcendental functions can come in quite 
different disguises. Might the left or the right hand side of (||) be a preferable 
representation? This question cannot be answered satisfyingly. A holonomic 
recurrence equation like (Q), defining the same transcendental function s(n), is 
probably the simplest way to describe a function of a discrete variable, since 
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it postulates how the values of the function can be calculated iteratively. Not 
only is this a quite efficient way to calculate the values of s(n), but moreover 
it is preferable to either of the two representations given in (^|), since it gives a 
unique representation scheme. This is what a normal form is about. 

As a further example, we consider the function (a, (3, 7 e No, z, M, d £ R + ) 

v(n r _ 1 ^+/3+ 7 r(a+/3+7-^)r(d/2- 7 )r(a + 7 -d/2)r(/3+7-d/2) 

1 ,P,1) [ 1 T{a)T(P)T{d/2)T{a + [i + 2 1 -d)M a +i 3 +~<- d 



■2FA 



a + (3 + 7 - d , a + j-d/2 
a + /3 + 2-f-d 



(2-F1 here represents Gaufi's hypergeometric function, see fill, C hapter 15), which 
plays a role for the computation of Feynman-diagrams |l5|El, for which Zeil- 
berger's algorithm generates the holonomic recurrence equation 

= {a + /3-d + j) (2a-d + 2j) V(a,/3,-y) 

+ aM(2a + 2/3-2d + 4j-2z-4az-2/3z + 3dz-4:jz) V(a + l,(3,j) 
+ 2 a (1 + a) M 2 {z - 1) z V{a + 2, /?, 7) 

and analogous recurrence equations with respect to the variables j3 and 7 (see 
[p4|). These, in particular, can be used for numerical purposes. 

Note that for the application of Zeilberger's algorithm our Mathematica 
program uses the Paule-Schorn implementation | 5jJ . For the current example, 
the output is given by 

In [18] : = HolonomicRE[(-l) " (alpha+beta+gamma) *Gamma[alpha+beta+gamma-d] * 
Gamma [d/2-gamma] *Gamma [alpha+gamma-d/2] *Gamma [beta+gamma-d/2] / 
(Gamma [alpha] *Gamma [beta] *Gamma [d/2] * 
Gamma [alpha+beta+2*gamma-d] *M~ (alpha+beta+gamma-d) ) * 
Hypergeomet r i c 2F 1 [alpha+bet a+gamma-d , alpha+gamma-d/ 2 , 
alpha+beta+2*gamma-d,z] , alpha ,V] 

Out [18]= (alpha + beta - d + gamma) (2 alpha - d + 2 gamma) V [alpha] + 

> alpha M (2 alpha + 2 beta - 2 d + 4 gamma - 2 z - 4 alpha z - 

> 2 beta z + 3 d z - 4 gamma z) V[l + alpha] + 

2 

> 2 alpha (1 + alpha) M (-1 + z) z V[2 + alpha] == 



^1 am indebted to Jochcm Fleischer who informed me about a misprint in formula (31) of 
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4 Holonomic Systems of Several Variables 



In |Q , Zeilberger considered the more general situation of functions F of several 
discrete and continuous variables. If we have d variables, and d (essentially in- 
dependent) mixed homogeneous linear (partial) difference-differential equations 
with polynomial coefficients in all variables are given for F, then F is called a 
holonomic system (compare H). In most cases these holonomic equations 
together with suitably many initial values declare F uniquely. 

In particular, we concentrate on the situation, when the given system of 
holonomic equations is separated, i.e. each of them is cither an ordinary differ- 
ential equation or a pure recurrence equation. These representing holonomic 
equations can be generated by the method described in § ^ whenever F is given 
in terms of sums and products of primitive functions. 

For example, the Legendre polynomials F(n,x) — P n (x) (|Q, Chapter 22) 
form a holonomic system by their holonomic differential equation 

{x 2 - l)F"(n,x)+2xF' (n, x) — n(l + n)F(n, x) = (8) 

and their holonomic recurrence equation 

(n + 2)F(n + 2,x) - (3 + 2n)xF(n + 1, x) + {n+ l)F(n, x)=Q, (9) 

together with the initial values 

F(0,0) = 1, F(1,0) = 0, F'(0,0) = 0, F'(1,0) = 1. (10) 

Equations (p|)-(|l0|) therefore build a sufficient algebraic, even polynomial struc- 
ture to represent the functions P n (x) as we shall see now. 

If we interpret the (partial) differentiations and shifts that occur as opera- 
tors, and the representing system of holonomic equations as operator equations, 
then these form a polynomial equations system in a noncommutative polyno- 
mial ring. For a continuous variable x with differential operator D given by 
DF(n,x) — F'(n,x), the product rule implies D(xf) — xDf = /, and hence 
the commutator rule Dx — xD = 1 is valid. Similarly for a discrete variable n 
with the (forward) shift operator N given by NF(n,x) — F(n + l,x), we have 
N(nF(n, x)) - nNF(n, x) = (n + l)F(n + 1, x) - nF(n + l,x) = F(n + l,x) = 
NF(n, x), and therefore the commutator rule Nn — nN — N. Similar rules are 
valid for all variables involved, whereas all other commutators vanish. 

The transformation of a holonomic system given by mixed holonomic differ- 
ence-differential equations represents an elimination problem in the noncom- 
mutative polynomial ring considered, that can be solved by noncommutative 
Grobner basis methods (§, @, @, g|, @, ||-|§), ||). 

Hence, we need the concept of a Grobner basis. If one applies GauB's al- 
gorithm to a linear system, the variables are eliminated iteratively, resulting in 
an equivalent system which is simpler in the sense that it contains some equa- 
tions which are free of some variables involved. Note that connected with an 
application of GauB's algorithm is a certain order of the variables. 
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The Buchbcrgcr algorithm is an elimination process, given a certain term 
order for the variables (a variable order is no longer sufficient), with which a 
polynomial system (rather than a linear one) is transformed, resulting in an 
equivalent system (i.e., constituting the same ideal) for which the terms that 
are largest with respect to the term order, are eliminated as far as possible. Note 
that — in contrast to the linear case — the resulting equivalent system may con- 
tain more polynomials than the original one. Such a rewritten system is called 
a Grobner basis of the ideal generated by the polynomial system given. It turns 
out that Buchberger's algorithm can be extended to the noncommutative case 
that we consider here |l9|j as long as the rewrite process using the commutator 
does not increase the variable order. 

As an example, we consider F(n, fc) — (^) in which case we have the Pascal 
triangle relation F{n + l,k + 1) = F(n,k) + F(n,k + 1), together with the 
pure recurrence equation (n + 1 — k)F(n + l,k) — {n + l)F(n,k) — with 
respect to n, say. These equations read as (if iV — 1 — K)F(n, k) — 0, and 
((n + 1 — k)N — (n + l))F(n, k) = in operator notation, K denoting the shift 
operator with respect to k. Therefore we have the polynomial system 

KN-l-K and (n + 1 - k)N - (n + 1) . (11) 



The Grobner basis of the left ideal generated by ( |1 1| ) with respect to the lexi- 
cographical term order (fc, n, K, N) is given by 

{(jfe + l)if + k - n, (n + 1 - k)N - (n + l),KN - 1 - if } , 

i.e., the elimination process has generated the pure recurrence equation 

(fc + l)F(n, k + 1) + (fc - n)F(n, fc) = 

with respect to fc. 

We used the REDUCE implementation |2£| for the noncommutative 
Grobner calculations of this article, but I would like to mention that there is also 
a Maple package Mgfun written by Chyzak JToj (to be obtained from 
|http : //pauillac . inria. f r/ algo/libraries/libraries .html#Mgf un ) which 



can be used for this purpose. 

As another example, we consider the Legendre polynomials. In operator 
notation the holonomic equations (||)-(^|) constitute the polynomials 

(x 2 -l)D 2 + 2xD-n(l + n) and (n + 2)N 2 - (3 + 2n)xN + (n + 1) . (12) 

The Grobner basis of the left ideal generated by (|l2|) with respect to the lexi- 
cographical term order (D, N, n, x) is given by 

|(.t 2 - 1)D 2 + 2xD - n(l + n), 

(l+n)ND- (l + n)xD- {l + nf, (13) 
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(x 2 -l)ND-(l + n)xN+(l + n), (14) 
(1 + n){x 2 - 1)D - (1 + n) 2 N + x(l + nf , (15) 

(n + 2)7V 2 - (3 + 2n)x7V + (n + 1)} . 

After the calculation of the Grobner basis, for better readability I positioned 
the operators £> and ./V back to the right, so that the equations can be easily 
understood as operator equations, again. By the term order chosen, the Grobner 
basis contains those equations for which the D-powers are eliminated as far as 
possible, and (p^)-(fl5|) correspond to the relations 

P; i+1 (x)=xP; i (x) + (l + n)P n (x) , 

(x 2 - l)P; i+1 (x) = (1 + n) (xP n+1 (x) - P n (x)) , 

(x 2 - 1)1* {x) = (1 + n) {P n+ i{x) - xP n (x)) (16) 

between the Legendre polynomials and their derivatives. 

If we are interested in a relation between the Legendre polynomials and 
their derivatives that is a;-free (which is of importance for example for spectral 
approximation, see ||), we choose the term order (x,D,N,n) to eliminate x 
in the first place, and obtain a different Grobner basis containing the rc-free 
polynomial 

-{n + 2)(n + 1)D - (2n + 3)(n + 2)(n + 1)JV + (n + 2)(n + 1)N 2 D 

equivalent to the identity 

(2n+l)P n (x)=P^ +1 (x)-P^_ 1 (x) 

for the Legendre polynomials (see e.g. ||, formula (2.3.16)). 

Here, we present the REDUCE output for the above examples: 

1: load ncpoly; 

2 : nc_setup ({D,NN,n,x}, {NN*n-n*NN=NN , D*x-x*D=l} , left) ; 

3: pl: = (x~2-l)*D~2+2*x*D-n*(l+n)$ "/„ differential equation 

4: p2: = (n+2)*NN~2-(3+2*n)*x*NN+(n+l)$ recurrence equation 

5: nc_groebner ({pi ,p2}) ; 

2 2 2 2 
{d *x - d - 2*d*x - n - n, 

2 

d*nn*n - d*n*x - d*x - n - n, 
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2 

d*nn*x - d*nn - nn*n*x - 2*nn*x + n + 1 , 



2 2 2 2 

d*n*x - d*n + d*x - d - nn*n + n *x - x, 

2 

nn *n - 2*nn*n*x - nn*x + n + 1} 

6 : nc_setup ({x,D,NN,n}, {NN*n-n*NN=NN , D*x-x*D=l} , left) ; 

7: result :=nc_groebner({pl ,p2}) ; 

2 2 2 2 

result := {x *d + 2*x*d - d - n - n, 

2 2 
x*d*nn - d*nn *n + d*n + d + 2*nn*n + 2*nn*n + nn, 

2 

x*d*n + x*d - d*nn*n + n + 2*n + 1, 
2 

2*x*nn*n + x*nn - nn *n - n - 1 , 

2 2 2 2 3 2 

d*nn *n - d*nn *n - d*n - 3*d*n - 2*d - 2*nn*n - 3*nn*n - nn*n} 

8 : nc_setup ({n , x , NN , D} , {NN*n-n*NN=NN , D*x-x*D=l} , left) ; 

9 : nc_compact (part (result , 5) ) ; 

2 

- (2*n + 3)*(n + 2)*(n + l)*nn + (n + 2)*(n + l)*nn *d - (n + 2)*(n + l)*d 

We see, therefore, that by the given procedure new relations (between the bino- 
mial coefficients, and between the derivatives of the Legendre polynomials) can 
be discovered. The generation of derivative rules like (|l6|), and the algorithmic 
work with them is described in [E3| . 



5 Holonomic Sums and Integrals 

Analogously, with the method in the last section, holonomic recurrence equa- 
tions for holonomic sums can be generated. Note that the idea to use recurrence 
equations for the summand to deduce a recurrence equation for the sum is origi- 
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nally due to Sister Celine Fasenmyer (Q-Q, see Q, Chapter 14). Zeilberger 
[ [l2| brought this into a more general setting. 
Consider for example 

fc=0 fc=0 ^ ' 

then by the product algorithm, we find the holonomic recurrence equations 
(n-k+ l)F(n + 1, k) - (1 + n)F(n, k) = 

and 

(2 + k) 2 F(n,k + 2)-(3 + 2k)(n-k-l)xF(n,k + l) + {n-k)(n-k-l)F(n,k) = 

for the summand F(n, k). The Grobner basis of the left ideal generated by the 
corresponding polynomials 

(n-k+l)N-(l+n) and (2+k) 2 K 2 -(3+2k)(n-k-l)xK+(n-k){n-k-l) 

with respect to the lexicographical term order (fc, iV, n, K) contains the fc-free 
polynomial 

{2+n) 2 K 2 N 2 -K(2 + n)(3 + 2n)(K+x)N+(l+n)(2+n)(l + K 2 +2Kx) , (17) 

which corresponds to a fc-free recurrence equation for F(n, k). We use the order 
(fc, N,n, K) because then fc-powers are eliminated as far as possible (since we 
like to find a fc-free recurrence), and A^-powers come next in the elimination 
process (since the recurrence equation obtained should be of lowest possible 
order). 

Because all shifted sums 

s{n) = F(n, k)=^2 F(n, k + 1) = ^ F(n, k + 2) 

generate the same function s(n), and since summing the fc-free recurrence equa- 
tion is equivalent to setting K = 1 in the corresponding operator equation 
(check!), the substitution K = 1 in ( p"7| ) generates the valid holonomic recur- 
rence equation 

(2 + n)s(n + 2) - (3 + 2n)(l + x)s(n + 1) + 2(1 + n)(l + x)s(n) = 
for s(n). 

In the general case, we search for a fc-free recurrence equation contained in 
a Grobner basis of the corresponding left ideal with respect to a suitably chosen 
weighted |3(J (or lexicographical (fc, N,n, K)) term order. For example, the 
elimination problems described in |45| are automated by this procedure. 
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On the other hand, it turns out that in many cases the holonomic recurrence 
equation derived is not of the lowest order. In the next section, we will discuss 
how this problem can be resolved. 

Note that by a similar technique, holonomic integrals can be treated gj. To 
find a holonomic equation for 



I(y) := / F(y,x)dx 

a 

for holonomic F(y, x) with respect to the discrete or continuous variable y, cal- 
culate the Grobner basis of the left ideal constituted by the holonomic equations 
of F(y, x) with respect to a suitably chosen weighted or the lexicographical term 
order (x, D y ,y, D x ). We search for an x-iree holonomic equation £ contained 
in such a Grobner basis. In case, that F(y, a) — F(y 7 b) = 0, and enough 
derivatives of F(y, x) with respect to x vanish at x — a and x — b, by partial 
integration it follows that the holonomic equation valid for I(y) is given by the 
substitution D x = into £ (see (|]). 
As an example, we consider 



I(n) := / e x H n (x) dx , 



H n {x) denoting the Hermite polynomials. The method of § g yields the holo- 
nomic polynomials 

2 (1 + n) + N 2 - 2xN and D 2 + 2 (1 + n) + 2 x D 

for the integrand. Note that since H n (x) is an odd function for odd n, it is 
immediately clear that I(n) = in this case. However, what about even values 
of nl 

The Grobner basis of the corresponding left ideal contains the two x-free 
polynomials 

N 2 + ND and Nn + nD + D 

so that setting D = we get for I{n) the recurrence equation I(n + 1) = 0. 
Indeed, this proves that I(n) = for n > 1. 

As another example, we consider the Abramowitz functions ([0, 27.5)) 

oo 

A(„,:=/*»e--^. 



(I 



By the method in § |2| for the integrand F(n,y,x) = x n e x y l x we get the 
three holonomic polynomials 

x — N , — n x + x 2 D x + 2 x 3 — y and 1 + x D y . 
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Using the term order (x, D y , y, D x ), the differential equation 

y A'"(n, y)-(n- 1) A"(n, y) + 2A(n, y) = , 

and using (x, N, n, D), the recurrence equation 

2A(n + 3,y)- (n + 2) A(n + l,y) -yA(n,y) = 

is automatically generated by the given approach (compare Jy, (27.5.1), (27.5.3)). 
Finally, we mention that similarly an identity like ([EJ, (11.4.28)) 



-a x „m—l 



X m - L J n {bx)dx 



T (n/2 + m/2) b n 

2 n+l a n+m T (n + 1) 



n/2 + to/2 
n + 1 



4a 2 



(18) 

(l-Fi representing Kummer's confluent hypergeometric function) for the Bessel 
function is proved by the calculation of the common holonomic recurrence equa- 
tion 







- (n + 3) (n + m) 6 2 /(n) 
+2 (n + 2) (4 a 2 n 2 + 16 a 2 n + 12 a 2 - b 2 m + b 2 ) I(n + 2) 
+ (n + 1) (n + 4 - to) 6 2 /(n + 4) 

for the left and right hand sides of (|l8|). Note that Zeilberger's algorithm is not 
directly applicable to the right hand side, but the extended version of [^3| gives 
the result. 



6 Noncommutative Factorization and Holonomic 
Normal Form 

Note that neither the sum and product algorithms of § || nor Zeilberger's algo- 
rithm or its extension |^3f , nor the algorithms for holonomic sums and integrals 
of § H can guarantee to present the holonomic equation J\f of lowest order, and 
therefore the normal form searched for. 

In (^9|lI a Grobner basis based factorization algorithm was introduced for 
polynomials in noncommutative polynomial rings given by Lie bracket commu- 
tator rules. This method is implemented in ]2q] . Given an expression /, and 
a holonomic equation V of order to of /, one may find the normal form M of 
/ using this factorization algorithm by generating the right factors of the non- 
commutative polynomial p corresponding to V, and checking if any of them, Q, 
say, (having order I < m, say) and to — I derivatives (shifts) of Q are satisfied 
by / at a certain initial point. In the affirmative case, Q is compatible with /, 
and corresponds to a valid holonomic equation for /. 

3 Due to a severe bicycle accident of Herbert Melenk, this paper is still unfinished. 
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To present some examples, we consider Zeilberger's algorithm first. An ex- 
ample for which Zeilberger's algorithm does not generate the holonomic recur- 
rence equation of lowest order is given by the sum (see e.g. pl[) 




for which the holonomic equation 

2 (2n + 3) s n+2 + 3 (5n + 7) s n+1 + 9 (n + 1) s n = (19) 

is generated. Note that there is an algorithm due to Petkovsek |3^] to find 
all hypergeometric solutions of holonomic recurrence equations which could be 
used as next step. However, we may also proceed as follows: The corresponding 
noncommutative polynomial 2(2n + 3)N 2 + 3(5n + 7)N + 9(n + 1) is factorized 
by implementation |28| as 

2(2n + 3)N 2 + 3(5n + 7) N + 9(n + 1) = ((An + 6) N + 3(n + 1)) (N + 3) . 

The right factor N + 3 corresponds to the holonomic recurrence equation 

Sn+i + 35„ = , (20) 

which, together with the initial value Sq — sq = 1 uniquely defines a sequence 
(S n ) n £fi . Since S\ = —3 turns out to be compatible with the given sum 




and since ( |20|) implies ( p9| ) (right factor!), the sequence s n , which is the unique 
solution of (|19|) with so = 1 and s\ — —3, must equal S n . From (pp[), however, 
the closed form s n = (—3)" follows. 

Similarly, for any particular d e N, d > 3, the identity 

S^G) (?)-<-<>" 

can be established, for whose left hand side Zeilberger's algorithm generates a 
recurrence equation of order d — 1 (see |3lJ ) . 

Whereas Petkovsek's algorithm finds hypergeometric solutions of holonomic 
recurrence equations as in the example, and therefore not only verifies identities, 
but generates closed-form results, our approach is more general in the follow- 
ing sense. Factorizations with polynomial coefficients of ordinary holonomic 
differential equations (see [Q, J35| for other methods) as well as of any mixed 
holonomic difference-differential equation can be calculated. 
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We give an example of that type for the sum algorithm: Consider the differ- 
ence of successive Gegenbauer polynomials h(x) = C^T \ (x) — C„ (x) that 

were used in j2(|. Here the summand f(x) := C« 1 (x) satisfies the holonomic 
equation 

(x 2 ~l)f"(x) + (n-n 2 )f(x) = 0, 
and the sum algorithm yields the fourth order equation 

(x 2 - l) 2 h""(x)+4x (x 2 -l) h"'(x)-2 (n 2 -l){x 2 -l) h"(x)+n 2 (n 2 -l) h(x) = 

for h(x). The implementation ]2§| ] finds (besides others) the noncommutative 
factorization 

((a; 2 - 1) D 2 + (1 + x) D - n 2 ) ((x 2 - 1) D 2 - (1 + x) D + (1 - n 2 )) 

of the corresponding noncommutative polynomial, whose right factor 

(x 2 - 1) D 2 - (1 + x) D + (l- n 2 ) 

turns out to be compatible with the given function hix). That is, the corre- 
sponding differential equation and two derivatives thereof are satisfied by h(x), 
at x = 1. Therefore the holonomic normal form of h(x) is the corresponding 
differential equation 

(1 - n 2 ) h{x) - (1 + x) h'(x) + {x 2 - 1) h"(x) = 

that was a tool in (26). This result can also be obtained by the method given 
in 

To evaluate the integrals 

oo 

/„ := J x n er x2 H n {x)dx , 

— oc 

we may deduce the holonomic system 

iV 2 -2:zi 2 7V + 2(l + n)a; 2 

and 

x 2 D 2 + 2x(a; 2 - n) D + (n + n 2 + 2x 2 ) 

of the integrand. The Grobner basis of this system with respect to the weighted 
lexicographical order with weights (3, 1,0,0) for {x, N,n, D) (i.e. the term x is 
considered larger than TV 3 , whereas x is smaller than TV 4 , and any power of n 
and D is smaller than x and N) contains an x-free polynomial, which when 
evaluated at D = yields 

P(n,N) = (n + 5)(n + 4)(n + 3)N 3 - (3n + 7)(n + 5)(n + 4)(n + 3)N 2 

+(3n + 5)(n + 5)(n + 4)(n + 3)(n + 2)N (21) 
-(n + 5)(n + 4)(n + 3)(n + 2)(n + l) 2 



19 



corresponding to a recurrence equation of order three. 

On the other hand, P(n, N) obviously has the trivial (commutative) factor- 
ization 

P(n,N) = (n+5)(n+4)(n+3)(iV 3 -(3n+7)7V 2 +(3n+5)(n+2)7V-( n+ 2)(n+l) 2 ^) 
and the remaining right factor can be represented as 

iV 3 -(3n+7)7V 2 +(3n+5)(7i+2)iV-(n+2)(n+l) 2 = (N-n-2)(N-n-l)(N-n-l) 

(note that finds four different right factors). This leads to the valid recur- 
rence equation I n +i = {n + l)I n that together with the initial value 

oo 

Io = J eT x dx = y/n 

— oo 

gives finally I n = y/nnl. 
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